%one can choose the following parameter values

A=*;

v_theta = 1; %this is normalization
v_3_def = *;    
v_2_def = *;  
v_1_def = *;
kappa_def = *;  

v_z_def= *;
mu_def = *;

phi_def =0; % the default value of \phi is 0

%%%%%%%%%%%%%%%%%%%%%%%%%%%

% the parameters are from the following ranges, respectively.
% A \in [1.5, 2.5);
% v_3_def, v_2_def, and v_1_def \in [0.05, 0.15) 

% kappa_def\in [1.45, 1.55)

%v_z_def\in [0.05, 0.15)
%mu_def \in [0.15, 0.25)

%%%%%%%%%%%%%%%%%%%%%%%%%%%
v_3 = v_3_def;  
v_2=v_2_def;
v_1=v_1_def;

kappa=kappa_def; 
v_z = v_z_def; v_z3=v_z;
mu = mu_def;
phi=phi_def;

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures 1 and 2: S and S_annual return depending on v_z and \mu
% Figures 3 and 4: L and L_annual return depending on v_z and \mu

I = 20;
vz_min = 0.00001; vz_max = 0.2;
MVZ = vz_min:(vz_max-vz_min)/(I-1):vz_max;

J=20;
mu_min = 0.0001; mu_max = 0.5;
MU  = mu_min:(mu_max-mu_min)/(J-1):mu_max;

S = zeros(I, J);  S_annual	=zeros(I, J); COV_E= zeros(I, J);
L = zeros(I, J);  L_annual 	= zeros(I, J);  

i=1;
while i<I+0.5	
	v_z = MVZ(i);	v_z3=v_z; 
	j=1;
	while j<J+0.5		
		mu = MU(j); 
		

		
        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
		S(i, j)	= rev; S_annual(i, j)	= rev_annual;   
		COV_E(i, j) =cov12_23;
		
		A=A_L;
		[rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
		L(i, j) = mom;  L_annual(i, j) = mom_annual;  		 
		 
		j=j+1;
	end
	i=i+1;	
end


mu = mu_def;
v_z = v_z_def; v_z3=v_z; 

Y = MVZ; X = MU; 


figure(1); Z = -S_annual;
mesh(X, Y, Z);

zticks([-0.1, -0.05, 0, 0.05, 0.1, 0.15, 0.2 0.25]); 

%zticklabels({'-10%', '-5%' '0', '5%', '10%', '15%', '20%',  '25%'});

   a=[cellstr(num2str(get(gca,'ztick')'*100))]; % Create a vector of '%' signs
   pct = char(ones(size(a,1),1)*'%');           % Append the '%' signs after the percentage values
   new_zticks = [char(a),pct];                  % 'Reflect the changes on the plot
   set(gca,'zticklabel',new_zticks)
   
xlabel('$\mu$','Interpreter','Latex');
ylabel('$\nu_z$','Interpreter','Latex');
title('$\displaystyle  -{ {\mathcal{S}}}/{|Y_{{\mathcal{S}}}|}$','Interpreter','Latex');
xlim([0, max(X)]);
ylim([min(Y), max(Y)]);
saveas(gcf,'fig1_S_annual','epsc');

figure(2); Z = L_annual;
mesh(X, Y, Z);

   a=[cellstr(num2str(get(gca,'ztick')'*100))]; % Create a vector of '%' signs
   pct = char(ones(size(a,1),1)*'%'); % Append the '%' signs after the percentage values
   new_zticks = [char(a),pct];% 'Reflect the changes on the plot
   set(gca,'zticklabel',new_zticks)

xlabel('$\mu$','Interpreter','Latex');
ylabel('$\nu_z$','Interpreter','Latex');
title('$\displaystyle   {\mathcal{L}}/{|Y_{ {\mathcal{L}}}|}$','Interpreter','Latex');
xlim([0, max(X)]);
ylim([min(Y), max(Y)]);
saveas(gcf,'fig2_L_annual','epsc');


figure(4); 
Z = COV_E;
mesh(X, Y, Z);
xlabel('$\mu$','Interpreter','Latex');
ylabel('$\nu_z$','Interpreter','Latex');
title('$  {\rm Cov}_E$','Interpreter','Latex');
xlim([0, max(X)]);
ylim([min(Y), max(Y)]); 
saveas(gcf,'fig4_covE','epsc');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This is for the default value for the S_annual


[rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
rev_annual



A=A_L;
[rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
mom_annual
diff= mom_skipping - mom

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figure 5: S, S(2), and S(3) depending on v_z, conditional on \mu=0 and \mu=0.3

mu=0.2;
I = 40;
v_z_min =  0.00001; v_z_max =  0.2;
X = v_z_min:(v_z_max-v_z_min)/(I-1):v_z_max;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

S = X; S2 = X; S3 = X;
REF_S1 = X; REF_S2 = X; REF_S3 = X;  

i=1;
while i<I+0.5	
	v_z = X(i); v_z3=v_z;	

        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual, ref_s1, ref_s2, ref_s3] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
        S(i) =rev; 
        S2(i)=(i1+i2)/2;     
        S3(i) =i3;  
        
        REF_S1(i) =ref_s1; 
        REF_S2(i)=ref_s2;     
        REF_S3(i) =ref_s3;  
	i=i+1;	
end

figure(3); 
plot(X, S, '-', X, S2, '+', X, S3, 'x');
ylim([-0.07, 0.05] );
xlabel('$\nu_z$','Interpreter','Latex');
leg1=legend(['$\mathcal{S}$' ], ['$\mathcal{S}_{(2)}$' ], ['$\mathcal{S}_{(3)}$' ],  'Location', 'southwest' );
set(leg1,'Interpreter','latex');
%title('Panel B: $\mu=0.2$','Interpreter','Latex');
saveas(gcf,'fig3_S_longlag','epsc');


mu=0.2;
I = 40;
v_z_min =  0.00001; v_z_max =  0.4;
X = v_z_min:(v_z_max-v_z_min)/(I-1):v_z_max;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
REF_S1 = X; REF_S2 = X; REF_S3 = X;  

i=1;
while i<I+0.5	
	v_z = X(i); v_z3=v_z;	

        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual, ref_s1, ref_s2, ref_s3] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);

        REF_S1(i) =ref_s1; 
        REF_S2(i)=ref_s2;     
        REF_S3(i) =ref_s3;  
	i=i+1;	
end

figure(11); 
plot(X, REF_S1, '-', X, REF_S2, '+', X, REF_S3, 'x');
%ylim([0, 0.05] );
xlabel('$\nu_z$','Interpreter','Latex');
leg1=legend(['${\rm Cov}(P_1-P_0, E_1(P_2-P_1))$' ], ['${\rm Cov}(P_1-P_0, E_1(P_3-P_2))$' ], ['${\rm Cov}(P_1-P_0, E_1(P_4-P_3))$'], 'Location', 'best' );
set(leg1,'Interpreter','latex');
%title('Panel B: $\mu=0.2$','Interpreter','Latex');
saveas(gcf,'ia_fig11_different_S_longlag','epsc');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Figure 7: S, Cov_E, and Cov_2E depending on v_z, conditional on \mu=0 and \mu=0.3

mu=0.2;
I = 40;
v_z_min =  0.00001; v_z_max =  0.2;
X = v_z_min:(v_z_max-v_z_min)/(I-1):v_z_max;

COV_E = X; COV_2E = X; Y=X;
mu = 0.2;  
i=1;
while i<I+0.5	
	v_z = X(i);		v_z3=v_z; 


        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
        S(i)=rev;

        COV_E(i) =cov12_23; COV_2E(i)=cov12_34;
	i=i+1;	
end

figure(5);
plot(X, S, X, COV_E, '+', X, COV_2E, 'x');
%ylim([-0.04, 0.04] );
xlabel('$\nu_z$','Interpreter','Latex');
%title('Panel B: $\mu = 0.2$','Interpreter','Latex');
leg1=legend(['$\mathcal{S}$' ], ['${\rm Cov}_E$' ], ['${\rm Cov}_{2E}$' ], 'Location', 'southwest');
set(leg1,'Interpreter','latex');
saveas(gcf,'fig5_covE_longlag','epsc');


mu=mu_def; v_z = v_z_def; v_z3=v_z;  


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 8: on v_z3

%%%%%%%%%%%%%%%%%%%%%%%%%%%

I = 40;
vz_min = 0.05; %0.000001; %0.07;
vz_max = 0.15; %5;
MZ3 = vz_min:(vz_max-vz_min)/(I-1):vz_max;

S = zeros(1, I); L = zeros(1, I); L_skipping = zeros(1, I);
C2334 = zeros(1,I);
AAA =zeros(1, I); BBB= zeros(1, I);
i=1;
while i<I+0.5	
	v_z3 = MZ3(i);	

	
    [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
	S(i)	= rev;

	A=A_L;
    
     [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual, ref_s1, ref_s2, ref_s3, cov01_12, cov23_34, aaa, bbb] = main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
	L(i) = mom;
	L_skipping(i)	= mom_skipping;  

    C2334(i)=cov23_34;
	
    AAA(i)=aaa;
    BBB(i)=bbb;

	i=i+1;	
end

v_z = v_z_def; v_z3=v_z;  



figure(6);
X = MZ3 +mu^2*v_z; 

subplot(3,1,1); 
Y = S;
plot(X, Y);
xlabel('$\nu_{z_3}+\mu^2 \nu_{z_2}$','Interpreter','Latex');
title('Panel A: $ {\mathcal{S}}$','Interpreter','Latex');
xlim([min(X), max(X)]);

subplot(3,1,2);
Y = L;
plot(X, Y);
xlabel('$\nu_{z_3}+\mu^2 \nu_{z_2}$','Interpreter','Latex');
title('Panel B: $ {\mathcal{L}}$','Interpreter','Latex');
xlim([min(X), max(X)]);

subplot(3,1,3);
Y = L_skipping;
plot(X, Y);
xlabel('$\nu_{z_3}+\mu^2 \nu_{z_2}$','Interpreter','Latex');
title('Panel C: $ {\mathcal{L}^*}$','Interpreter','Latex');
xlim([min(X), max(X)]);
saveas(gcf,'fig6_LS_on_vz3','epsc');

figure(17);
Y=C2334;
plot(X, Y);
xlabel('$\nu_{z_3}$','Interpreter','Latex');
title('$ {\rm Cov}(P_3-P_2, P_4-P_3)$','Interpreter','Latex');
xlim([min(X), max(X)]);
saveas(gcf,'fig6_cov2334_on_vz3','epsc');

figure(18);
subplot(2,1,1);
Y=AAA;
plot(X, Y);
xlabel('$\nu_{z_3}$','Interpreter','Latex');
title('$dP_2/dF_2$','Interpreter','Latex');
xlim([min(X), max(X)]);

subplot(2,1,2);
Y=BBB;
plot(X, Y);
xlabel('$\nu_{z_3}$','Interpreter','Latex');
title('$ dP_2/dz_2$','Interpreter','Latex');
xlim([min(X), max(X)]);
saveas(gcf,'fig6_P2_F2z2_on_vz3','epsc');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures  9, 10: on \phi

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is to plot Ss against mu2

I=40; x_min = 0.0001; x_max = 0.3;
X  = x_min:(x_max-x_min)/(I-1):x_max;
mu =0.3;
S = X; S2 = X; S3=X; L=X;
i=1;
while i<I+0.5	
	phi = X(i); 	


    [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
	S(i)	= rev; S2(i)=(i1+i2)/2; S3(i) =i3;
    
    A=A_L;
    [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
    L(i) = mom;
   
	i=i+1;	
end

figure(7); 
plot(X, L, '-');
xlabel('$\phi$','Interpreter','Latex');
title('$ {\mathcal{L}}$','Interpreter','Latex');
%leg1=legend(['$\mathcal{L}$'] );
%set(leg1,'Interpreter','latex',  'Location', 'best');
saveas(gcf,'fig7_L_on_phi','epsc');

figure(8); 
plot(X, S, '-', X, S2, '+', X, S3, 'x');
xlabel('$\phi$','Interpreter','Latex');
leg1=legend(['$\mathcal{S}$' ], ['$\mathcal{S}_{(2)}$' ], ['$\mathcal{S}_{(3)}$' ] );
set(leg1,'Interpreter','latex',  'Location', 'best');
saveas(gcf,'fig8_S_on_phi','epsc');

mu=mu_def; phi=phi_def;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 11: L and S depending on kappa

I = 40;
kappa_min =  1.5; kappa_max =  2.5;
X = kappa_min:(kappa_max-kappa_min)/(I-1):kappa_max;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

S = X; L=X;  
i=1;
while i<I+0.5	
	kappa = X(i);	


        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);

       S(i) =rev;
		
		A=A_L;
		[rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
		L(i) = mom; 		 
       
	i=i+1;	
end

kappa=kappa_def; 

figure(11); 
yyaxis left
plot(X,S, 'b-')
xlabel('$\kappa$','Interpreter','Latex'); % X-axis label for the first plot
ylabel('${\mathcal{S}}$','Interpreter','Latex'); % Y-axis label for the first plot

ax = gca;

% Set the color of the x-axis and y-axis to black
ax.XColor = 'k'; % Black color for x-axis
ax.YColor = 'b'; 

yyaxis right
plot(X,L, 'r--')
ylabel('${\mathcal{L}}$','Interpreter','Latex'); % Y-axis label for the first plot

ax = gca;

% Set the color of the x-axis and y-axis to black
ax.YColor = 'r'; % Black color for y-axis

leg1=legend(['$\mathcal{S}$' ], ['$\mathcal{L}$' ] );
set(leg1,'Interpreter','latex',  'Location', 'best');

saveas(gcf,'fig_LS_on_kappa','epsc');


figure(31); 

plot(X,L)
xlabel('$\kappa$','Interpreter','Latex'); % X-axis label for the first plot
title('${\mathcal{L}}$','Interpreter','Latex'); % Y-axis label for the first plot

saveas(gcf,'fig_L_on_kappa','epsc');
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 6: L and S depending on v_2

I = 40;
v_2_min =  0.4; v_2_max = 0.6;
X = v_2_min:(v_2_max-v_2_min)/(I-1):v_2_max;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

S = X; L=X;  COV_E=X;
i=1;
while i<I+0.5	
	v_2 = X(i);	

        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);

       S(i) =rev;
       		COV_E(i) =cov12_23;
		
		A=A_L;
		[rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
		L(i) = mom; 	
   
	i=i+1;	
end

v_2=v_2_def;
figure(12); 
yyaxis left
plot(X,S, 'b-')
xlabel('$\nu_2$','Interpreter','Latex'); % X-axis label for the first plot
ylabel('${\mathcal{S}}$','Interpreter','Latex'); % Y-axis label for the first plot

ax = gca;

% Set the color of the x-axis and y-axis to black
ax.XColor = 'k'; % Black color for x-axis
ax.YColor = 'b'; 

yyaxis right
plot(X,L, 'r--')
ylabel('${\mathcal{L}}$','Interpreter','Latex'); % Y-axis label for the first plot
leg1=legend(['$\mathcal{S}$' ], ['$\mathcal{L}$' ] );
set(leg1,'Interpreter','latex',  'Location', 'best');
ax = gca;
% Set the color of the x-axis and y-axis to black
ax.YColor = 'r'; % Black color for y-axis
saveas(gcf,'fig_LS_on_v2','epsc');


figure(32); 
plot(X,L);
xlabel('$\nu_2$','Interpreter','Latex'); % X-axis label for the first plot
title('${\mathcal{L}}$','Interpreter','Latex'); % Y-axis label for the first plot
 
saveas(gcf,'fig_L_on_v2','epsc');




figure(13);
plot (X, COV_E);
xlabel('$\nu_2$','Interpreter','Latex');
title('$ {\rm Cov}_E$','Interpreter','Latex');
xlim([min(X), max(X)]);
saveas(gcf,'fig_covE_on_v2','epsc');

figure(14);
plot (X, S);
xlabel('$\nu_2$','Interpreter','Latex');
title('$\mathcal{S}$','Interpreter','Latex');
xlim([min(X), max(X)]);
saveas(gcf,'fig_S_on_v2','epsc');


figure(15);
plot(X, COV_E-S);
xlabel('$\nu_2$','Interpreter','Latex');
title('$  {\rm Cov}_E - \mathcal{S}$','Interpreter','Latex');
xlim([min(X), max(X)]);
saveas(gcf,'fig_covE_deduct_S_on_v2','epsc');

figure(16);
plot (X, COV_E, '-', X, S, '--');
xlabel('$\nu_2$','Interpreter','Latex');
leg1=legend(['${\rm Cov}_E$' ], ['$\mathcal{S}$' ] );
set(leg1,'Interpreter','latex',  'Location', 'best');
xlim([min(X), max(X)]);

saveas(gcf,'fig_covE_and_S_v2','epsc');




%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure cov_E on v_2 and v_3

I = 20;
v2_min = 0.4; v2_max = 0.6;
MV2 = v2_min:(v2_max-v2_min)/(I-1):v2_max;

J=20;
v3_min = 0.4; v3_max = 0.6;
MV3  = v3_min:(v3_max-v3_min)/(J-1):v3_max;

S = zeros(I, J);  S_annual	=zeros(I, J); COV_E= zeros(I, J);
L = zeros(I, J);  L_annual 	= zeros(I, J);  

i=1;
while i<I+0.5	
	v_2 = MV2(i);	
	j=1;
	while j<J+0.5		
		v_3 = MV3(j); 
		

		
        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
		
		COV_E(i, j) =cov12_23;
		 
		 
		j=j+1;
	end
	i=i+1;	
end

v_2 = v_2_def;
v_3 = v_3_def;

Y = MV2; X = MV3; 

figure(19); 
Z = COV_E;
mesh(X, Y, Z);
xlabel('$\nu_3$','Interpreter','Latex');
ylabel('$\nu_2$','Interpreter','Latex');
title('$  {\rm Cov}_E$','Interpreter','Latex');
xlim([min(X), max(X)]);
ylim([min(Y), max(Y)]); 

saveas(gcf,'fig_covE_v2_v3','epsc');


%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure : 01-12 depending on v_1

I = 40;
v_1_min =  1; v_1_max = 10;
X = v_1_min:(v_1_max-v_1_min)/(I-1):v_1_max;
S = X;  COV_E=X;

i=1;
while i<I+0.5	
	v_1 = X(i);	

        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual, ref_s1, ref_s2, ref_s3, cov01_12, cov23_34] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);

       S(i) =rev;
       COV_E(i) =cov01_12;
	i=i+1;	
end
v_1=v_1_def;

figure(21); 

yyaxis left
plot(X,COV_E, 'b-')
xlabel('$\nu_1$','Interpreter','Latex'); % X-axis label for the first plot
ylabel('$ {\rm Cov}(P_1-P_0,P_2-P_1 )$','Interpreter','Latex'); % Y-axis label for the first plot
xlim([min(X), max(X)]);
ax = gca;

% Set the color of the x-axis and y-axis to black
ax.XColor = 'k'; % Black color for x-axis
ax.YColor = 'b'; 

yyaxis right
plot(X,COV_E-S, 'r--')
ylabel('$ {\rm Cov}(P_1-P_0,P_2-P_1 )-\mathcal{S}$','Interpreter','Latex'); % Y-axis label for the first plot
xlim([min(X), max(X)]);
leg1=legend('${\rm Cov}(P_1-P_0,P_2-P_1 )$' , '${\rm Cov}(P_1-P_0,P_2-P_1 ) - \mathcal{S}$'  );
set(leg1,'Interpreter','latex',  'Location', 'best');
ax = gca;
% Set the color of the x-axis and y-axis to black
ax.YColor = 'r'; % Black color for y-axis

saveas(gcf,'fig_covE_0112_deduct_S_on_v1','epsc');



%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure : 12-23 depending on v_2

I = 100;
v_2_min =  0.4; v_2_max = 10;
X = v_2_min:(v_2_max-v_2_min)/(I-1):v_2_max;
S = X;  COV_E=X;

i=1;
while i<I+0.5	
	v_2 = X(i);	

        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual, ref_s1, ref_s2, ref_s3, cov01_12, cov23_34] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);

       S(i) =rev;
       COV_E(i) =cov12_23;
	i=i+1;	
end
v_2=v_2_def;

figure(22); 

yyaxis left
plot(X,COV_E, 'b-')
xlim([min(X), max(X)]);
xlabel('$\nu_2$','Interpreter','Latex'); % X-axis label for the first plot
ylabel('$ {\rm Cov}(P_2-P_1,P_3-P_2 )$','Interpreter','Latex'); % Y-axis label for the first plot

ax = gca;

% Set the color of the x-axis and y-axis to black
ax.XColor = 'k'; % Black color for x-axis
ax.YColor = 'b'; 

yyaxis right
plot(X,COV_E-S, 'r--')
xlim([min(X), max(X)]);
ylabel('$ {\rm Cov}(P_2-P_1,P_3-P_2 )-\mathcal{S}$','Interpreter','Latex'); % Y-axis label for the first plot
leg1=legend('${\rm Cov}(P_2-P_1,P_3-P_2 )$' , '${\rm Cov}(P_2-P_1,P_3-P_2 ) - \mathcal{S}$'  );
set(leg1,'Interpreter','latex',  'Location', 'best');
ax = gca;
% Set the color of the x-axis and y-axis to black
ax.YColor = 'r'; % Black color for y-axis

saveas(gcf,'fig_covE_1223_deduct_S_on_v2','epsc');



%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure : 23-34 depending on v_3

I = 100;
v_3_min =  0.4; v_3_max = 10;
X = v_3_min:(v_3_max-v_3_min)/(I-1):v_3_max;
S = X;  COV_E=X;

i=1;
while i<I+0.5	
	v_3 = X(i);	

        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual, ref_s1, ref_s2, ref_s3, cov01_12, cov23_34] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);

       S(i) =rev;
       COV_E(i) =cov23_34;
	i=i+1;	
end
v_3=v_3_def;

figure(23); 

yyaxis left
plot(X,COV_E, 'b-')
xlim([min(X), max(X)]);
xlabel('$\nu_3$','Interpreter','Latex'); % X-axis label for the first plot
ylabel('$ {\rm Cov}(P_3-P_2, P_4-P_3)$','Interpreter','Latex'); % Y-axis label for the first plot

ax = gca;

% Set the color of the x-axis and y-axis to black
ax.XColor = 'k'; % Black color for x-axis
ax.YColor = 'b'; 

yyaxis right
plot(X,COV_E-S, 'r--')
xlim([min(X), max(X)]);
ylabel('$ {\rm Cov}(P_3-P_2, P_4-P_3)-\mathcal{S}$','Interpreter','Latex'); % Y-axis label for the first plot
leg1=legend('${\rm Cov}(P_3-P_2, P_4-P_3)$' , '${\rm Cov}(P_3-P_2, P_4-P_3) - \mathcal{S}$'  );
set(leg1,'Interpreter','latex',  'Location', 'best');
ax = gca;
% Set the color of the x-axis and y-axis to black
ax.YColor = 'r'; % Black color for y-axis

saveas(gcf,'fig_covE_2334_deduct_S_on_v3','epsc');



%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure: L and S3 depending on kappa

I = 100;
kappa_min =  1.5; kappa_max =  2.5;
X = kappa_min:(kappa_max-kappa_min)/(I-1):kappa_max;
%%%%%%%%%%%%%%%%%%%%%%%%%%%

S3 = X; L=X;  
i=1;
while i<I+0.5	
	kappa = X(i);	


        [rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);

       S3(i) =i3;
		
		A=A_L;
		[rev, i1, i2, i3, rev_annual, mom, mom_skipping, cov12_23, cov12_34,  mom_annual] =main_func(A, v_theta, v_3, v_2, v_1, kappa, v_z, v_z3, mu, phi);
		L(i) = mom; 		 
       
	i=i+1;	
end

kappa=kappa_def; 

figure(24); 
yyaxis left
plot(X,S3, 'b-')
xlabel('$\kappa$','Interpreter','Latex'); % X-axis label for the first plot
ylabel('${\mathcal{S}}_{(3)}$','Interpreter','Latex'); % Y-axis label for the first plot

ax = gca;

% Set the color of the x-axis and y-axis to black
ax.XColor = 'k'; % Black color for x-axis
ax.YColor = 'b'; 

yyaxis right
plot(X,L, 'r--')
ylabel('${\mathcal{L}}$','Interpreter','Latex'); % Y-axis label for the first plot

ax = gca;

% Set the color of the x-axis and y-axis to black
ax.YColor = 'r'; % Black color for y-axis

leg1=legend(['$\mathcal{S}_{(3)}$' ], ['$\mathcal{L}$' ] );
set(leg1,'Interpreter','latex',  'Location', 'best');

saveas(gcf,'fig_L_S3_on_kappa','epsc');

